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Abstract 



We investigate the one-dimensional quantum XXZ model in the presence 



S3 . 

of diagonal disorder. Recently the model has been analyzed with the help 

O 

of field-theoretical renormalization group methods, and a phase diagram has 
been predicted. We study the model with exact diagonalization techniques 
up to chain lengths of 16 sites. Using finite-size scaling methods we estimate 
critical exponents and the phase diagram and find reasonably good agree- 
ment with the field-theoretical results, namely, that any amount of disorder 
destroys the superfluidity for XXZ anisotropy A between —1/2 and 1, while 
the superfluidity persists to finite disorder strength for —1 < A < —1/2 and 

then undergoes a Kosterlitz-Thouless type transition. 
PACS numbers: 64.70.Pf, 05.30.Jp, 05.70.Jk, 75.10.Nr 
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I. INTRODUCTION 



The two focal points of many-body physics in recent times are the study of strongly 
interacting systems and the exotic phenomena induced by the presence of disorder. The 
repeated refinements of the analytic and numerical tools and the accumulation of experience 
in these major areas have set the stage for the quest of understanding strongly interacting 
disordered systems. For fermionic problems the attack started in the field of disordered 
semiconductors jT],3. On the other hand, bosonic studies were delayed by the haunting 
puzzles raised by their classical counterpart, the spin-glass problem; and by the relative 
inaccessibility of such systems experimentally. The past few years brought breakthroughs 
in both of these directions. A fairly comprehensive picture of the spin-glass state has been 
developed ||, and elaborate techniques established the dirty superconductors f|-§], helium 
in vycor 0H, and quantum spin chains as well controllable experimental realizations of 
the model, posing many challenges for theoretical studies. 

We concentrate on the phenomena at zero temperature, as the most profound differences 
between the ordered and disordered systems manifest themselves at that point. At T = 
the ordered models undergo a quantum phase transition, which occurs as a parameter of the 
Hamiltonian is tuned across some critical value. In this case quantum fluctuations drive the 
transition instead of the usual thermal ones. Ordered d dimensional quantum systems are 
equivalent to a corresponding d + 1 dimensional classical systems and as such their critical 
phenomena are well understood. Whereas for disordered quantum systems the classical 
analogues are much less worked out - e.g., the McCoy- Wu model flO| is the single exactly 
solved model - and as such their study demands genuinely new theoretical approaches. The 
subject of our paper is the numerical study of this competition between quantum ordering 
tendencies and the disruptive effects of disorder on the example of the one-dimensional 
disordered XX Z model. 

The ordered one-dimensional XX Z model with spin 1/2, described by the Hamiltonian: 
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H = J2 (SfSf +1 + SVS? +1 + A SiSfa) (1) 

8=1 

is one of the most-studied quantum systems and many of its properties are well known 



lT|.[T2[. For the values of the anisotropy parameter — 1 < A < 1 the system develops quasi- 
long-range order (quasi-LRO). Expectation values of the spin-operators vanish, but the 
spin-spin correlation functions decay only as a power-law in the ground state. For A > 1 
the excitation spectrum is gapped and long range antiferromagnetic order of the Ising-type 
is established, and similarly for A < — 1 ferromagnetic Ising long range order sets in. 

As it is well known, the spin 1/2 XX Z model is equivalent to a lattice gas of hard 
core bosons, which is thought of as an approximate representation of 4 He []13| , ^4| . The z- 
component of total magnetization and boson particle number are related by M z + L/2 = N^. 
The above antiferromagnetic phase corresponds to the solid phase of helium. Ordering in 
the XY plane maps onto the superfluid phase. 

The disorder will be represented by adding a random magnetic field to the Hamiltonian: 

L 

Hrandom 2 ^ tl^ . (2) 
i=l 

where {hi ) = 0, {{hi) 2 ) = D. In the equivalent hard-core boson problem this corresponds to 
the inclusion of a random site-energy. For other types of disorder, such as random bonds, see 



Refs. |T5| and |16| . The effect of disorder in the Ising regimes is well understood. According 
to the Imry-Ma argument the long range order is destroyed by the addition of a weak random 
magnetic field in the ^-direction, however it is expected to persist in the presence of a weak 



random exchange term |T7 



The quasi-LRO regime has been first studied by Giamarchi and Shulz (GS) in the pres- 
ence of disorder |18| , who developed a powerful scaling scheme for the problem. They utilize 
the Haldane representation for the bosons ||19|| , when writing down the effective action: 

S = Jdx dr^[(d T ®) 2 + (d x <S>) 2 } + [Ti(x)d x * + p(x)e i$ + p*{x)e~ i % (3) 

Here $ is the phase field, representing the bosons and rj and p are the Fourier components 
of the disorder fields at momenta k pa and k pa 7rpo, respectively, where po is the average 



density of the bosons. In the general picture of disordered systems, the backscattering (i.e., 
the large momentum component) is driving the localization phenomena and this expectation 
is borne out by explicit calculation in the present case as well. In Eq. (0) k is a spin-stiffness 
of the ordered system. It can be related to the original parameters by analyzing the Bethe- 
ansatz solution of the problem for zero disorder to arrive at 



k = - -cos" 1 A). (4) 
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Upon integrating out the disorder one arrives at an action similar to that of the sine-Gordon 
problem: 
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S eS = J2j dx dr-[(d T ^ a ) 2 + (d x $ c 

n , 

-D ^2 / dx dr dr'cos[^ a (x, r) - $p(x, t')], (5) 

a,/3=l 

where the a and (3 sums are over the n replicas. The infrared singularities are then taken 
care of by a renormalization group analysis, which yields: 

dD/dl = (3- \/2<kk)D, 

dn/dl = l -D, (6) 

where / = In b, and b is the scale change ratio. For small values of the randomness and 
interactions not too strong, so that k < k c = 1/67T, D renormalizes towards a line of fixed 
points at D — 0, and thus the quasi-LRO persists in the presence of disorder. However above 
that critical value of k c the disorder becomes a relevant operator, destroying the ordering 
tendencies already for arbitrary small values of D. The phase transition is analogous to 
the Kosterlitz-Thouless (KT) type. This scaling analysis can be viewed as the quantum- 
generalization of the Harris criterion. 



Recently Doty and D. Fisher gave an exhaustive study of the phase diagram ||15|| . By uti- 
lizing several scaling arguments in different parameter regimes, they constructed a schematic 
phase diagram for the case of weak random z— fields. In particular, the GS phase transition 



occurs at A = —1/2. For — 1/2 < A the disorder is relevant and in very small amounts it 
destroys the quasi-LRO. Nagaosa has also come to this conclusion p^| . This result is not 
surprising for the A = (XY-model) case which maps exactly onto noninteracting fermions 
in a disordered potential |£J and is well-known to become localized with infinitesimal disor- 
der. In the —1 < A < —1/2 regime the quasi-LRO is argued by Doty and Fisher to be stable 
against the disorder up to a finite value of D = D c . The disordered phase was identified as 



a "bose-glass" by M. Fisher et al. p4|. The main physical feature of the bose-glass phase is 
that all of its low-lying excitations are localized. Consequently, it has a vanishing superfluid 
density and a finite compressibility. In other words, the localization is achieved not by the 
opening of a gap in the spectrum, i.e. via the Mott scenario, but rather by the Anderson 
mechanism, which localizes the particles by the interference of their wavef unctions. In the 
localized state the spin-spin correlations decay exponentially with spatial separation. By in- 
tegrating the recursion relations Eqs. (|^) away from the Kosterlitz-Thouless critical regime 
the following relation is obtained for the correlation length in the region —1/2 < A < 1: 

Z~D-*% (7) 

where 

0, = (3- 1/21XK)- 1 (8) 

is a crossover exponent, k is that of the pure system (Eq. (fD) and so depends only on 
the anisotropy A. Eq. (§) is essentially the the scaling dimension of the Born-scattering 
amplitude. On the other hand, in the —1 < A < —1/2 region spin-correlations should still 
decay as a power-law for small D > 0, however the exponent of the spin-spin correlation 
function is modified by the presence of disorder. The KT transition occurs when disorder 
increases the stiffness to k — l/6n. In our work we set out to perform an extensive numerical 
survey of the above ideas. 
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II. NUMERICAL METHOD 



We use the exact diagonalization method to compute the properties of the disordered 
spin 1/2 XX Z model. The random fields h\ in Eq. are drawn from an independent 
and uniform distribution at each site, such that (/if) = 0, {{hi) 2 ) = D. For each realization 
of disorder the ground state \%Ijq) is found via an accelerated power method |2"5fl . Once the 
vector |^o) is computed, expectation values of observables, such as the ground state energy 
E and the spin-spin correlation functions, are determined. The chemical potential fi and 
compressibility K are found by taking discrete derivatives of the ground state energy with 
respect to the boson number N^: [i = OE/dN^ and K~ l = Ld 2 E/dN^ 2 . All results reported 
here are for the half-filled density {i.e. M z = 0) case: = L/2. 

The superfluid density - which is related to the helicity modulus - is computed from the 
formula |26| 



1 d 2 E 

Ps = 0) 

via finite differencing with respect to 6, where 9 is the angle of a phase twist applied at the 
boundary. A boson hopping to the right through the boundaries acquires a phase e l9 , while 
one hopping in the opposite direction acquires e~ l6 '. One may think of p s as a measure of 
the "degree of sensitivity to boundary conditions": in the quasi-LRO superfluid state phase 
coherence is long-ranged enough to yield a finite p s , whereas in the localized phase p s drops 
off exponentially with system size L. 

For each (A, D) pair, 300 to 5, 000 realizations of disorder are used for averaging. The 
system sizes we study are L = 4, 6, . . . , 14, 16, with the smallest number of realizations for 
the larger systems. 



III. ANALYSIS OF DATA 

First we checked the accuracy of our numerical procedure on the clean system. We 
calculated the spin-spin correlation function = {S^S* + SfS^) which should behave as 
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where 77 = 2hk. For a finite-size system one would expect the correlation at separation 
L/2 should go as Z/ _?? , and so one may estimate 77 by: 

ln[r(L72)/r(L/2)] 
' ~ ln[L'/L] 1 ' 



where L and V are two different system sizes [p7| . Indeed, using this method |23] we find 
an exponent within one percent of the value given by the Bethe-ansatz solution for nearly 
all values of anisotropy A (except near A = 1, the isotropic HAF point, where logarith- 
mic corrections complicate our extrapolation scheme). We also find these estimates for 77 



agree closely with that predicted from the thermodynamical quantities via 27T77 = 1/ y/p s K. 
Results for 77, p s , and the compressibility K are shown in Fig. (1). Note that while p s is 
nearly constant in the entire — 1 < A < 1 range, K varies strongly, and diverges at the 
isotropic ferromagnetic point A = — 1 (where all bosons occupation number sectors have 
the same ground state energy, thereby making the system infinitely compressible). The 
prediction [p||15| as to whether the system's quasi-LRO will be stable or unstable with 



respect to the addition of infinitesimal disorder depends only on the pure system quantity 
k = l/(47r 2 v / p s K). As the compressibility increases the superfluid phase correlations be- 
come stronger (i.e. decay more slowly) until the point is reached where a D c > is required 
to drive p s to zero. 

We started the study of the disordered system by computing the universal scaling function 
for the superfluid density. Utilizing the relation between the current-current correlation 
function and p s f29|] it is straightforward to derive the finite-size scaling ansatz: 



Ps (L,0~C~ {d+z) Ps(L/0, (11) 



where p s (x) is a universal function. As shown by Giamarchi and Shulz ||18|| , the dynamical 
critical exponent z is one in one dimension, thus in fact the superfluid density itself is 
expected to be a universal function in our case 

We determined p s for roughly 1000 points |HJ in the (L, A, D) parameter space. Typical 
data for p s is shown in the insets to Figs. (2a), (2b), (2c), and (2d), and later in Figs. 



(8a) and (8b). Before delving into the detailed analysis, some qualitative statements can be 
made. In the large disorder regime (i.e. yD > 0.30), p s vanishes quickly with L, whereas 
for D — it clearly extrapolates to a non-zero value, indicating the presence of a phase- 
transition. Thus the system is much more susceptible to the addition of disorder than is 
its two-dimensional (2D) analog. The 2D case with A = and Hamiltonian given by the 
sum of Eqs. ([!]) and @ was considered by Runge f25|, where it was found that the critical 



value of disorder was \fTT c = 1.3. We will see below that for the one dimensional (ID) case 
in the whole range — 1 < A < 1, \fTT c is never bigger than 0.1 or so WA - This means for 
— 1 < A < 1 the value of the D parameter in the effective action of Eq. (|5]) that destroys 
quasi-LRO is no larger than 0.04. Such small values of critical disorder, however, are not 
unexpected: 1) the superffuid phase is only power-law correlated as opposed to 2D which 
possesses true LRO, and 2) D c is expected to be zero for A < — 1 and for —1/2 < A, so it 
is not too surprising that D c cannot become large in the remaining interval. 
Power-law regime: — 1/2 < A < 1 

To our knowledge the first numerical work on the spin 1/2 XX Z model in the presence 
of disorder is due to Nagaosa who studied < A and utilized a transfer matrix method 
based on the Suzuki- Trotter breakup to deduce the finite temperature properties of a long 
(L = 200 and 10 Trotter time slices) chain. Nagaosa predicted scaling relations and obtained 
very good scaling functions for the superfluid and charge density wave susceptibilities as 
1/kT = (3 —>■ oo. Our work complements his, as we explicitly have (3 = oo and study 
finite system sizes. We focus on the quantities p s ,K, and £. For some properties it is 
advantageous to have (3 = oo since this is where the quantum critical phenomenon is more 
naturally described, that is to say, the phase transition occurs in the ground state as a 
parameter in the Hamiltonian is varied. 

1) For — 1/2 < A < 1 the correlation length is predicted to depend on the disorder according 
to Eqs. (0), that is, £ ~ D~^ s . Therefore as a first method, motivated by this form and by 
Eq. flTTD, we consider p s (L,D) as a function of x = (D — D c ) L 1 ^ 3 , and then choose the 
parameters D c and <p s that collapse the different p s curves best. This is done via a non-linear 
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(x 2 ) least squares fit. Since the data points almost never have the same value of x as defined 
above, spline fits are used to interpolate data points to common x values. In most cases data 
from L = 8 or 10 to 16 is used since the smaller system sizes tend to not be in the scaling 
regime. In Figs. (2a, 2b, 2c, and 2d) examples of the data collapse for A = 0, 0.5, —0.25, 1.0 



are presented |33fl. The data evidently scale well even out to large D where p s becomes small. 



For these values of A, the optimum value of the fitting parameter D c was zero ||34j| . In Fig. 
(3) we show best fit values of <p s as circles for the different anisotropics A. For comparison 
we display Doty and Fisher's prediction for the exponent <p s . The data appears consistent 
with the results of the scaling analysis Hl8|Jl5| , namely, that the critical value of the disorder 
is zero and that the crossover exponent is given by Eq. (]8|). The accurate extraction of <p s 
and D c is hampered, however, as the KT region at A = —1/2 is approached. This is not 
too surprising as we are using small systems to discern a crossover between different forms 
of critical behavior rather than the critical behavior alone. To underpin the above results 
we conducted several further analyses, discussed in the following. 

2) An unconstrained extraction of the correlation length can be performed from the super- 
fluid density data. We plot p s against ln(L) for the different values of the disorder D. Then 
we shift each curve horizontally to maximize overlap with its 'neighboring' curve possessing 
the nearest value of D. This process is repeated for all .D-curves to obtain a sequence of 
shifts, an example is shown in Fig. (4a). As p s is presumed to be a function exclusively 
of L/£, the amount of shift to maximize the overlap of the curves with Di and D 2 is given 
by hi(j;(Di) / £(D 2 )) . From all the shifts one can build up, to within an overall constant, 
the function £(D). This method was used by Yoshida and Okamoto in their study of the 
ID spin 1/2 XX Z model in the presence of an alternating bond perturbation (which has 
a remarkable degree of similarity to the present disordered system) [p5| . Once the optimal 
shifts are performed the scaling function p s (x) is obtained again. 

Our results for A = and 0.5 are shown in Figs. (4a) and (4b), and demonstrate again 
that the data scale nicely. The inset to Fig. (4a) shows the original unshifted data, p s vs. 
ln(L). We show the extracted correlation length £(D) for different A in Fig (5). Note the 
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potential power of this method to yield information on the bulk correlation length even when 
it is much larger then the system sizes considered. A log-log plot of £(D) against D displays 
a linear region from which one can extract an estimate of </> s ; the results are displayed as 
crosses in Fig. (3). The values of <p s are close to, and, unfortunately, not more accurate 
than those computed with method 1) of constrained data collapse. This is not surprising 
since the two analyses are closely related. 

We note in passing that we have also extracted a "system size dependent" correlation 
length from an equation similar to Eq. (|T0|) : 



ffl * >™W2)] (12) 

where V is a system size close to L. We have found £l(-D) does behave as D~^ B for a range 
of D, however, when £i becomes comparable to L it, of course, deviates from the D~^ s form. 
The finite-size scaling ansatz for £l(D) should be ^l{D)/L = F(£(D)/L), where F(x) is a 
universal function and £(D) is the bulk correlation length. Hence the shifting method of 
the previous paragraph should also be applicable to this data as well. For other possible 
definitions of £l(D) see Ref. P 



3) In our third method we consider p s cLS cL function of d = yD for the different values of 
system size L. By the symmetry /if — > —hf, p s must be an "even" function of d (note that 
d rather than D enters directly into the Hamiltonian). In particular, dp s /dd = at d = 
and at d — oo. Therefore, dp s /dd takes on a maximum that should be in the transition 
region since for all A we expect a jump in p s across d c for the infinite system. Let d*(L) 
be the point of maximum slope of p s (L,d). As L — > oo, d*(L) should tend to the location 
of the bulk critical point, d c = \[D~ C - Finite size scaling, Eq. (PD, indicates this maximum 
should occur at d*(L) - d c ~ if d c ^ 0, and d*(L) ~ 1/L 1 / 2 ^ if d c = 0. Therefore, 

we compute dp s (L, d)/dd, find the value d*(L) where it is maximal, and then plot d*(L) vs. 
L -i/2^_ pig _ ( 6 ) showg guch plotg for A = _ .38, -0.25, 0.0, 0.5, 0.9, 1.0. As can be seen in 

the figure, all the data is consistent with the transition occurring at D = (as d*(L) > 
the A = —0.38 curve must evidently bend upward), that is to say, disorder is a relevant 
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perturbation. 

4) Finally, we can proceed by accepting that D c = and attempt to extract <fr s for each A 
studied. Consider now p' s (0) = dp s (L, D)/dD at D = 0. This is the initial rate of change of 
the superfluid density with respect to the perturbation D. According to Eqs. (|7|) and ( |TT|) 
this object should behave as ~ L 1 /^. This method was introduced by Yoshida and Okomoto 
p5| . Note that in general the unlimited growth of p' s {0) with L is a further indication that 
D c = H36f . The the slope of p' s (0) vs L on a log-log plot should reveal the value of the 



4> s exponent. Our results are shown in Fig (7) |37|. The straight lines in the plots is the 
prediction of Doty and Fisher |i~5| . The <p s are extracted via a linear least squares fit and 
are plotted as squares in Fig. (3). There is very good agreement between our numerical 
results and the predictions of Ref. |L5|: they agree within 4%, 3%, and 2% at A = 0.0,0.5 



and 1.0, respectively. When A is close to —1/2, however, the analysis is less accurate since 
the critical phenomenon is evidently crossing over to a different form. On the basis of the 
above four different methods we can declare with good deal of certainty that the predictions 
of the scaling theory are confirmed by our data in the — 1/2 < A < 1 regime. 
Kosterlitz-Thouless regime: — 1 < A < — 1/2 

Next we study the —1 < A < —1/2 region, where it is predicted [0,0 that weak 
disorder is not relevant so that there is a finite region in the parameter space where the 
quasi-long-range-order survives. Here one extracts the dependence of the correlation length 
on the disorder by fully integrating the renormalization group Eqs. In this case we 

have to represent that the disorder D itself is strongly modified by the scaling. Close to the 



critical point the integration gives a Kosterlitz-Thouless-type formula [15 



i{D) ~ exp(A/y/D - D c ) (13) 

for D > D c . The Kosterlitz-Thouless (KT) behavior is difficult to extract even in clean 
systems, and even more so for disordered ones. Thus we confine ourselves to show that 
the data are consistent with a KT form, and attempt to estimate the phase boundary. We 
analyze the critical behavior with the same methods as above. 
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1) We first construct the universal scaling function for the superfluid density. The only 
difference is that £ is now given by Eq. fll3"|) , so we plot p s against Lexp(— A/yJD — D c ). 
and fit A and D c to minimize \ 2 . The scaling functions found for A = —0.8 and —0.6 are 
shown in Figs. (8a) and (8b). Good scaling is evidently achieved, however, the value of D c 
is poorly determined. Roughly, \fTT c is found to be about 0.1 or smaller, with an error bar of 
the same magnitude. We believe this problem is due to the very small value of D c , since the 
finite-size rounding region is actually a good deal larger than D c itself. This behavior can 
be contrasted with the 2D classical XY model on the square lattice. For that model, the 
"relatively large" Kosterlitz-Thouless T c can be located quite accurately from small systems. 
We have performed Monte Carlo simulations of the 2D classical XY model helicity modulus 
(p s ) to test our methods of locating Kosterlitz-Thouless phase transitions. Using systems 
with linear dimension up to L = 16 and the same finite-size scaling methods applied here 
(and below) we could readily locate T c to within 4% of its accurately established value |38 



It is unfortunate that the present model requires L much larger than 16 for accurate results 
39| , since this is near the limit of the exact diagonalization method. 



We try two fitting forms for £(-D) in addition to Eq. ([13]): (1) £ ~ exp(y4/v / c? — d c ), where 
d = v^TJ, and (2) £ ~ \D — D c \~ u . For D c ^ and sufficiently near the transition form 
(1) is equivalent to Eq (0). We observe for A = —0.8 and —0.6 the use of form (1) lowers 
the optimal \ 2 by about 30% from that we achieve from use of Eq. ([13]). This observation 
points to the fact that D c is so small, higher order terms in £ play a large role in describing 
the data. Use of form (2), £ <~ \D — D c \~ u , yields over a factor of two increase in the value 
of x 2 relative to that we find with Eq. (pT3|). This is somewhat promising because it at least 
suggests the mechanism for the transition in — 1 < A < —1/2 is of the Kosterlitz-Thouless 
form. Similarly, in the region < A < 1 we find the use of the KT £(-D) (Eq. (|T3"|)) fitting 
form yields \ 2 values 3 to 8 times larger than that from the power-law form for £(-D), thereby 
bolstering the belief that the power-law is the correct form for — 1/2 < A < 1. 

2) For completeness, we apply method 2) of the previous section involving p s vs InL 
shifting to obtain the scaling function to the —1 < A < —1/2 data. The scaling function 
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for A = —0.5 is shown in Fig. (9). Once again, although a reasonable scaling function is 
obtained, the extraction of D c values from the resulting £(D) proves difficult. The correlation 
length £(D) does appear to grow much more rapidly than the data for < A, which at least 
hints at the expected KT behavior. 
3) We apply the maximum slope method as well. Setting £ ~ L yields the finite-size 
scaling prediction d*(L) — d c ~ (1/lnL) 2 . Fig. (10) displays the d*(L) data for A = 
—0.95, —0.8, —0.6, —0.5. All extrapolate to values of y/~D~ c < 0.15, consistent with the values 
of D c crudely estimated with method 1) above. The value nearest the isotropic ferromagnetic 
point, A = —0.95, is the smallest with \fTT c ~ 0.03 and is consistent with the prediction of 



Doty and Fisher that D c — > as A — > — 1 from above [40|. However, this method yields 
its largest value of D c at A = — 1/2 which is inconsistent with the renormalization group 
prediction that D c = at the point where the pure system k = l/6n. To investigate this 
effect further, we apply the d*(L) vs. (1/lnL) 2 method to the remainder of the data (i.e. 
— 1/2 < A < 1 where the power-law form for £ and D c = are expected). The estimated 
\flT c values are plotted as squares in Fig. (11). The error bars for these points are estimated 
roughly by looking at the linear least squares fitting error for linear and for quadratic fits 
to d* vs (1/lnL) 2 , and also by extrapolating D*(L) = (d*(L)) 2 with similar forms. Near 
A = this method is consistent again with D c = as was shown previously using the 
power-law extrapolation. There is, in fact, enough curvature in the d*(L) vs. (1/lnL) 2 
plots for A = —0.25 and 0.0 to suggest that the simple extrapolations we perform here are 
inadequate. As mentioned before, in the predicted crossover regime at A = —1/2 severe 
finite-size errors should be expected. This is more than likely reflected in our overestimate 
of the size of the superfluid region with method 3). It could be that the almost linear curve 
viewed over a limited range of (1/lnL) 2 actually bends over for larger L. So, it is possible 
that the present extrapolation method for sufficiently large systems would yield D c = for 
A = — 1/2. At present it is impossible to know for certain whether or not this will happen. 
4) As a fourth method we develop a very different approach that makes explicit use of the 
renormalization group results. From our p s and K (compressibility) data for finite systems in 

13 



the presence of disorder we can compute the quantity k = l/4ii 2 ^p s K. The scaling analysis 
suggests that the bose-glass transition takes place when k = I/671", (i.e. when the power-law 
decay exponent rj = 1/3). However in the A < —1/2 regime the disorder renormalizes the 
value of ft, thus the pure system Bethe-ansatz-based Eq. (|]) cannot be applied. 

This method involves finding the point D*(L) where k(L, D*) = 1/67T, and extrapolating 
D*(L) to L = 00. Finite-size scaling indicates the correction should again go as (l/lnL) 2 . 
This k extrapolation method has, of course, D c = for —1/2 < A guaranteed (for large 
enough L). Since the field theoretic arguments [p^ , ^5|| are compelling [41], this method 



may prove to be the most accurate one near A = —1/2. The results of this extrapolation 
method are shown as diamonds in Fig. (11). Error estimates are performed as above in 
3) by examining a number of extrapolation forms along with the statistical error. It is 
promising that at A = —0.8 this method gives the same value of D c as did the maximum 
slope method 3). However, at A = —0.95 the ^-extrapolation method predicts a somewhat 
larger value of D c than method 3), although they evidently agree within large error bars. 
Near A = — 1 there is a large renormalization of the compressibility K as disorder is turned 
on. For example, at A = —0.95 the pure system K is near 9, whereas in the transition 
region it is reduced to around 4. This rapid variation introduces a large extrapolation of 
d*(L) to the thermodynamic limit. As an aside, we mention that all of our numerical data 
strongly imply a finite compressibility K at the transition points in accordance with the 
general predictions of the transition to the Bose Glass phase [^J[ . 



5) The final method is very similar to the previous one. It utilizes the interesting property 
of the KT recursion relations that the finite-size corrections to k in Eqs. (^) o,t the transition 
point are universal [0. This result may be derived by expanding k near the transition point 



clS Ki — 1/6-7T + e, then Eqs. (^) become: 

dD/dl = aeD 
de/dl = l -D, (14) 
where a = 18tt. These yield dD/de = 2ae which may be integrated to give D = ae 2 + const. 
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The critical manifold has const = 0, which upon integrating from a starting length scale Iq 
to the system size I = In L gives 

^)-^ = f"-« < 15) 

or e(l) = —2/(al) as I gets large. Thus at the transition point D c we expect the finite-size 
formula 

k(L,D c ) = ^-—\— + •••. (16) 

Note that the presence of the cutoff / = In L is (necessarily) asymptotically independent of 
the fraction of L one selects, since ln(aL) = InL + In a ~ InL. This method has been used 
for the 2D classical XY model to locate T c ||38|| , and so we attempt to use it here on our 
disordered quantum XXZ model. Fig. (12) shows the 2tt[k(L, D) — l/Gn] vs l/ln(L) data 
for A = —0.8 along with the straight line prediction of Eq. ([16]). A quadratic curve has been 
added to suggest the possible asymptotic behavior (i.e. critical manifold) of the D = 0.125 
curve. It is reassuring that this value of D c is consistent with and close to the values obtained 
in methods 3) and 4). These three methods show that our results for A = —0.8 are indeed 
strongly suggestive of the Kosterlitz-Thouless scenario. 

As a final indication that the superfluid phase persists to finite disorder D, consider 
Fig. (13), where we plot the log of the spin-spin correlation function (S*Sj + SfSj) versus 
ln(L) in the top panel and versus L in the bottom panel. The data is for A = —0.8 and 
\[T) = 0.100 and 0.225. In the top panel a power-law decay should yield a straight line, 
and, indeed, the \[~D = 0.100 data does this, whereas the \H5 = 0.225 curve drops off more 
quickly, indicating exponential decay. In the bottom panel the log-linear plot does suggest 
that the \fT) = 0.225 plot is approaching straight line behavior (i.e. exponential decay), 
whereas the \f~D = 0.100 still has significant upward curvature. These data are consistent 
with the previous p s , K analyses that indicate D c is around 0.125 or so. 

Before closing we mention that recently a Quantum Monte Carlo study of the bose- 
Hubbard model has been performed ||42|| , which is expected to exhibit the same critical 
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behavior and generic phase diagram as the present model. The existence of superfluidity in 
the presence of disorder was well established in that study. Indeed, in the soft core model of 
Ref. JEJ] (i.e. with coulomb repulsion U < oo), in general one expects a larger compressibility 
K than for the corresponding hard core (U = oo) model. Through rj oc l/y/p s K, this leads 
to smaller values of 77 and hence a wider range of stability of the superfluid phase. Since 
the model of Ref. |^2j corresponds to A = (i.e. no nearest neighbor interaction), there 
should exist a certain value of U at which the 77 of the pure system will be 1/3, and above 
this value of U infinitesimal disorder destroys the superfluidity. 

In conclusion we studied the quantum spin 1/2 XX Z model in one-dimension with 
diagonal disorder via exact diagonalization techniques. By employing different finite-size 
scaling methods we mapped out the phase diagram shown in Fig. (11). We found that 
weak disorder is relevant in the — 1/2 < A < 1 regime, or equivalently the critical value of 
the disorder is zero. Our estimate of the power-law exponent S was found in agreement 
with the predictions of Doty and Fisher |Lj]. For — 1 < A < — 1/2 the results suggest 
that a small, but finite disorder is needed to destroy the quasi-LRO. Thus at small disorder 
the superfluidity prevails as shown in Fig. (11). Our data in the transition region can be 
described by the Kosterlitz-Thouless form. Therefore, the overall picture emerging from our 
analysis is in agreement with the field theoretical and renormalization group predictions of 
Refs. H, m, and Q22|. 
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IV. FIGURE CAPTIONS 



Fig. (1) Pure system superfluid density p s , power-law decay exponent r], and compress- 
ibility K for the ID spin 1/2 XX Z model with z-anisotropy A. A = 1 corresponds to the 
isotropic Heisenberg antiferromagnet (HAF). Curves are the exact Bethe- ansatz solutions for 
the infinite system. Symbols denote the extrapolation of our L < 16 exact diagonalization 
results. p s and K were obtained from numerical finite differences, while 77 was extracted 
from the spin-spin correlation function (S^Sj + SfSj). 

Fig. (2) Finite size scaling function for the superfluid density as a function of the com- 
bination (D — D c )L l ^ s (that is simply related to £/L) for the optimal values of the fitting 
parameters D c and <p s (method 1). Insets show the original unsealed p s data. Anisotropies 
are: (a) A = 0.0 (XY model), (b) A = 0.5, (c) A = -0.25, and (d) A = 1.0 (HAF). 

Fig. (3) Exponent for the power-law divergence of the correlation length: £(D) ~ \D — 
D c \'^ a . Circles denote values determined by least squares fitting to D c and 4> s to best 
collapse the p s data (method 1). Crosses were obtain from the unconstrained £ determination 
(method 2) and have been displaced horizontally by 0.05 for clarity. Squares are extracted 
from the divergence of dp s /dD at D = as L — > 00 (method 4). The curve is the field 
theoretic renormalization group prediction of Ref. |15| (Eq. ([8])). 

Fig. (4) p s scaling function from unconstrained correlation length extraction method (2) 
[ 35| . Inset shows p s vs InL data before the shift. Straight lines connect the data in the inset 
with common disorder D values. Anisotropies are: (a) A = 0.0, (b) A = 0.5. 

Fig. (5) Correlation lengths £(D) computed via method 2 (p s vs InL shift method). Note 
that the extracted £ is much bigger than the maximum system size (L = 16). The A values 
are listed in the figure, where one notes that for fixed D, £ is a monotonically decreasing 
function of anisotropy A. 

Fig. (6) Finite-size extrapolation of the position of maximum slope d*(L) of p s (L, d) using 
the value of exponent cj) s given in Eq. @ (method 3). Extrapolation of d* to zero implies 
infinitesimal disorder destroys the superfluidity, i.e. disorder is a relevant perturbation. 
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Fig. (7) Log-log plot of the initial slope (with respect to D) of p s vs L for a number of 
anisotropies A (method 4). Straight lines have slope 1/0 S (with <p s from Eq. (|8|)) and are 
the predictions of the field theoretical treatment (Ref. ||15|| ). 

Fig. (8) p s scaling function based on the assumption of the Kosterlitz-Thouless (KT) 
correlation length divergence (Eq. (|TBD), in the region of A where the KT form is expected 
(method 1). Insets show original unsealed data as in Fig. (2). Anisotropies are: (a) 
A = -0.8, (b) A = -0.6. 

Fig. (9) A = —0.5 scaling function for p s from the unconstrained £(-D) determination 
(method 2), as in Fig. (4). 

Fig. (10) Extrapolation to L = oo of the position of maximum slope d*(L) of p s {L, d), for 
anisotropies A = —0.5, —0.6, —0.8, and —0.95 in the region where the Kosterlitz-Thouless 
transition is expected (method 3). 

Fig. (11) Disorder (D) - Anisotropy (A) phase diagram for the random field, one- 
dimensional, spin 1/2 XX Z model. The system is superfluid with quasi-LRO for the pure 
system D = with — 1 < A < 1. For -1/2 < A < 1 it is predicted |jl||T3| that infinites- 
imal disorder destroys the quasi-LRO. For — 1<A<— 1/2 it is predicted that there is a 
Kosterlitz-Thouless transition at finite D c > 0. The crosses (x) are our results from method 
(1) of constrained p s data collapse and are consistent with D c = 0; the squares (□) are 
from the maximum slope method (3) assuming the KT form of the transition; the diamonds 



(o) are from method (4) involving the stiffness criterion (T8|,[T5| k(D c ) = 1/6ty. The region 
" SF" indicating the superfluid phase in the presence of finite disorder is a semi-quantitative 
estimate of the phase boundary. It is a parabola constrained to vanish at A = — 1 and —1/2 
and go through our estimate at A = —0.8 [f5J. 

Fig. (12) Examination of the "flow" of the superfluid stiffness k = 1/ (Air 2 ^/ p s K) with 
increasing system size L (method 5). k c = l/fwr. The straight line is the renormalization 
group prediction of the finite-size correction at D c (Eq. (fCB])). The curved line meeting the 
data point for D = 0.125 is drawn to suggest the critical manifold: above this line the system 
flows to a localized state (k — > oo, p s — > 0), whereas below it the system flows to a quasi- 
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LRO/superfluid state with p s > and power-law decay exponent (rj < 1/3) renormalized 
by disorder. 

Fig. (13) Log-log and log-linear plots of (SfS* + SfSj) at separation L/2 for A = —0.8, 
and \/~D = 0.100 and 0.225. Top panel straight line behavior suggests the \/~D = 0.100 
point is in the quasi-LRO superfluid phase, while straight line behavior in the bottom panel 
suggests the \f~D = 0.225 point has exponentially decaying correlations, i.e. it is in a 
localized phase. 
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